// Figure 1

use "data/figure1/simulation_main.dta", clear

// Preliminary operations
label variable t "From 1 to 100"
egen t_demmaj = group(t) if dm==1
egen t_repmaj = group(t) if dm==0

sum t_repmaj
sum t_demmaj


// Sample with different fraction of dem years
g sample10 = 1 if t_demmaj<=5 | t_repmaj<=45
label var sample10 "10% dem years"

g sample20 = 1 if t_demmaj<=10 | t_repmaj<=40
label var sample20 "20% dem years"

g sample30 = 1 if t_demmaj<=15 | t_repmaj<=35
label var sample30 "30% dem years"

g sample40 = 1 if t_demmaj<=20 | t_repmaj<=30
label var sample40 "40% dem years"

g sample50 = 1 if t_demmaj<=25 | t_repmaj<=25
label var sample50 "50% dem years"

g sample60 = 1 if t_demmaj<=30 & t_demmaj!=.  | t_repmaj<=20 & t_repmaj!=.
label var sample60 "60% dem years"

g sample70 = 1 if t_demmaj<=35 & t_demmaj!=. | t_repmaj<=15 & t_repmaj!=.
label var sample70 "70% dem years"

g sample80 = 1 if t_demmaj<=40 & t_demmaj!=. | t_repmaj<=10 & t_repmaj!=.
label var sample80 "80% dem years"

g sample90 = 1 if t_demmaj<=45  & t_demmaj!=. | t_repmaj<=5 & t_repmaj!=.
label var sample90 "90% dem years"

foreach n in 10 20 30 40 50 60 70 80 90 {
	display "sample`n'"
tab dm if sample`n'==1
corr d m  if sample`n'==1
}


// Time dummies
quietly: tab t, g(dt)

// save intermediate dataset
save "data/figure1/simulation_main_intermediate.dta", replace


******************************************
*********ANALYSIS***********
*******************************************


scalar drop _all

// Settings
global kernel uniform
global y y_full


foreach sample in 10 20 30 40 50 60 70 80 90 {
use "data/figure1/simulation_main_intermediate.dta", clear

keep if sample`sample'==1

// select bandwidth
 rdbwselect $y dgap, kernel($kernel) 
	forval i = 5(1)20 {
	scalar h`i'=`i'*e(h_mserd)/10
	g bw`i'=0
	replace bw`i'=1 if dgap>-`=h`i'' & dgap<`=h`i''
	} 
	
	
	forval i = 5(1)20 {
		reg $y d dgap rdgap if bw`i'==1, rob 
			scalar bA`i'= _b[d]
			scalar seA`i'= _se[d]
		reg $y d m dt* dgap rdgap if bw`i'==1, rob 
			scalar bD`i'= _b[d]
			scalar seD`i'= _se[d]
		
		matrix a`i' =[h`i',bA`i',seA`i',bD`i',seD`i']
		} 
	
	matrix a=[a5 \	a6 \	a7 \	a8 \	a9 \	a10 \	a11 \	a12 \	a13 \	a14 \	a15 \	a16 \	a17 \	a18 \	a19 \	a20 ]
		
	svmat a
	rename a1 h
	label variable h "Bandwidth"
	
	rename a2 bA
	rename a3 seA
	g upA = bA+1.96*seA
	g lowA =  bA-1.96*seA
	
	rename a4 bD
	rename a5 seD
	g upD = bD+1.96*seD
	g lowD =  bD-1.96*seD
	

// plot effects as a function of bandwidth	
	line  bA upA lowA bD upD lowD h, xline(`=h10')  title(`: variable label sample`sample'', size(medium) color(black)) graphregion(color(white))   ///
			lcolor (black black black red red red) ///
			lstyle(solid shortdash shortdash solid shortdash shortdash solid shortdash shortdash) ///
			ylabel(0(0.1)0.6) ///
			legend(order(1 4 7) label(1 "Model A")  label(4 "Model B")  rows(1)  region(lcolor(white) lwidth(0)))
	

	graph save "results/sim_sample`sample'.gph" , replace
}
	
	scalar drop _all


 // combine all graphs
grc1leg results/sim_sample10.gph results/sim_sample20.gph results/sim_sample30.gph results/sim_sample40.gph results/sim_sample50.gph results/sim_sample60.gph results/sim_sample70.gph results/sim_sample80.gph results/sim_sample90.gph, legendfrom(results/sim_sample90.gph) graphregion(color(white)) 
graph display, ysize(18) xsize(14)

// save multi-graph
graph export "results/simulation.pdf", replace

// erase single graphs
foreach f in results/sim_sample10.gph results/sim_sample20.gph results/sim_sample30.gph results/sim_sample40.gph results/sim_sample50.gph results/sim_sample60.gph results/sim_sample70.gph results/sim_sample80.gph results/sim_sample90.gph {
	cap erase "`f'"
}

